#include "base.h"

int main(int argc,char **argv)
{

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@///

	  //////////////////
	  // Input   ///////
	  //////////////////

	  char kappa_sea[1024];
	  int conf_ini, conf_end, jump;
	  int Nconf, Nex, Nclust;
	  int Nmass,Nr;
	  int Nf;
	  int T,L;
          double beta, Uplaq;
	  float theta[4];
	  char dir_input[1024];
	  char dir_output[1024];
	  char gauge_fixing[1024];
	  char gauge_action[1024];

	

	  ////// Parsing command line

    	  for(int i=1; i<argc; i++) {
	    if(strncmp(argv[i],"-kappa_sea",10)==0)  sscanf(argv[i],"-kappa_sea=%s",kappa_sea);
	    if(strncmp(argv[i],"-conf_ini",9)==0)  sscanf(argv[i],"-conf_ini=%d",&conf_ini);
	    if(strncmp(argv[i],"-jump",5)==0)  sscanf(argv[i],"-jump=%d",&jump);
	    if(strncmp(argv[i],"-Nclust",7)==0)  sscanf(argv[i],"-Nclust=%d",&Nclust);
	    if(strncmp(argv[i],"-Nex",4)==0)  sscanf(argv[i],"-Nex=%d",&Nex);
	    if(strncmp(argv[i],"-Nmass",6)==0)  sscanf(argv[i],"-Nmass=%d",&Nmass);
	    if(strncmp(argv[i],"-Nr",3)==0)  sscanf(argv[i],"-Nr=%d",&Nr);
            if(strncmp(argv[i],"-Nf",3)==0)  sscanf(argv[i],"-Nf=%d",&Nf);
            if(strncmp(argv[i],"-dir_input",10)==0)  sscanf(argv[i],"-dir_input=%s",dir_input);
            if(strncmp(argv[i],"-dir_output",11)==0)  sscanf(argv[i],"-dir_output=%s",dir_output);
            if(strncmp(argv[i],"-gauge_fixing",13)==0)  sscanf(argv[i],"-gauge_fixing=%s",gauge_fixing);
            if(strncmp(argv[i],"-gauge_action",13)==0)  sscanf(argv[i],"-gauge_action=%s",gauge_action);
	    if(strncmp(argv[i],"-T",2)==0)  sscanf(argv[i],"-T=%d",&T); 
            if(strncmp(argv[i],"-L",2)==0)  sscanf(argv[i],"-L=%d",&L); 
	    if(strncmp(argv[i],"-beta",5)==0)  sscanf(argv[i],"-beta=%lf",&beta);            
            if(strncmp(argv[i],"-Uplaq",6)==0)  sscanf(argv[i],"-Uplaq=%lf",&Uplaq);

    	  }

	  Nconf=Nclust*Nex;
	  conf_end=conf_ini+jump*(Nconf-1);
          printf("==============================\n");
          printf("conf_ini=%d, conf_end=%d jump=%d\nNconf=%d, Nclust=%d, Nex=%d\n",conf_ini, conf_end,jump,Nconf,Nclust,Nex);
          printf("Nmass=%d\n",Nmass);
          printf("L=%d T=%d\n",L,T);
          printf("Nr=%d \n",Nr);
          printf("dir_input=%s \n",dir_input);
          printf("dir_output=%s \n",dir_output);
          printf("==============================");
	

	  //// Momenta set


	   char input_momdat[1024];
	   sprintf(input_momdat,"mom.dat");
    	   FILE *momdat=fopen(input_momdat,"r");
	   int Nindexmom;

	   fscanf(momdat,"%f %f %f %f",
           &theta[0],&theta[1],
           &theta[2],&theta[3]
           );
    	   fscanf(momdat,"%i",&Nindexmom);

          int ***iP=(int***)malloc(sizeof(int**)*(Nindexmom));
          for(int imom=0;imom<Nindexmom;imom++) iP[imom]=(int**)malloc(sizeof(int*)*(4));
          for(int imom=0;imom<Nindexmom;imom++) for (int idir=0; idir<4; idir++) iP[imom][idir]=(int*)malloc(sizeof(int)*(2));

	   for(int inm=0; inm<Nindexmom; inm++)
	    {
	      fscanf(momdat,"%i %i %i %i %i %i %i %i",
        	   &iP[inm][0][0],&iP[inm][0][1],
	           &iP[inm][1][0],&iP[inm][1][1],
	           &iP[inm][2][0],&iP[inm][2][1],
        	   &iP[inm][3][0],&iP[inm][3][1]
          	 );
   	    }
  
    	  fclose(momdat);
	  // Default values for testing
	  //iP[0][0][0]=0;iP[0][0][1]=3;iP[0][1][0]=0;iP[0][1][1]=2;iP[0][2][0]=0;iP[0][2][1]=2;iP[0][3][0]=0;iP[0][3][1]=2;
	  //iP[1][0][0]=4;iP[1][0][1]=7;iP[1][1][0]=2;iP[1][1][1]=3;iP[1][2][0]=2;iP[1][2][1]=3;iP[1][3][0]=2;iP[1][3][1]=3;
	 


	  int Nmom_per_set[Nindexmom];
	  for(int imom=0;imom<Nindexmom;imom++){
			Nmom_per_set[imom]=1;
			for(int idir=0; idir<4; idir++) Nmom_per_set[imom]=Nmom_per_set[imom]*(iP[imom][idir][1]-iP[imom][idir][0]+1);
	     }
	   int Nmom_tot=0;
	   for (int imom=0; imom<Nindexmom;imom++) Nmom_tot=Nmom_tot+Nmom_per_set[imom];

	   int Nmom;
	   Nmom=Nmom_tot;	

          // Momentum variables definition

              double *P2=(double*)malloc(sizeof(double)*(Nmom));
              double *SinP2=(double*)malloc(sizeof(double)*(Nmom));
              double **P=(double**)malloc(sizeof(double*)*(Nmom));
              for (int imom=0 ; imom<Nmom; imom++) P[imom]=(double*)malloc(sizeof(double)*(4));
              double **SinP=(double**)malloc(sizeof(double*)*(Nmom));
              for (int imom=0 ; imom<Nmom; imom++) SinP[imom]=(double*)malloc(sizeof(double)*(4));
              double *SinP4=(double*)malloc(sizeof(double)*(Nmom));
	      double *nn=(double*)malloc(sizeof(double)*(Nmom));
	  
	   // CALL RUTINE MOMENTUM


	      Momentum(iP,theta,P2,SinP2,P,SinP,SinP4,Nindexmom,Nmom_per_set,T,L);
	      compute_nn(iP,Nindexmom,nn);

	  //////////////////////////////////
          // IO files
          ///////////////////////////////

          char s0ff_base_input[1024];
          char s0ft_input[1024];
          char s0ff_base_output[1024];
          char s0ff_output[1024];

	  char file[1024];

        
          sprintf(s0ff_base_input,"%s/s0ft/s0ft",dir_input);
          sprintf(s0ff_base_output,"%s/FF_ALLP/ff",dir_output);


//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

           /////////////////////////////
	   // Read Fourier transform
	   //////////////////////////////
	
          //def
	  colorcolorspinspin ****s0FT_r=malloc(sizeof(colorcolorspinspin***)*(Nconf));
          for (int iconf=0; iconf<Nconf; iconf++) s0FT_r[iconf]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
          for (int iconf=0; iconf<Nconf; iconf++)  for(int imass=0; imass<Nmass; imass++) s0FT_r[iconf][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom)); 
	  for (int iconf=0; iconf<Nconf; iconf++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT_r[iconf][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));
	
	  int conf_to_read;
	  FILE *fin=NULL;


           printf("====================\n");
           printf("READING...\n");

           for (int r=0; r<Nr; r++)
	   {
           for (int iconf=0; iconf<Nconf; iconf++)
 	   {
	           conf_to_read=conf_ini+iconf*jump;
                   //printf("  ... conf %d\n",conf_to_read);
                   sprintf(s0ft_input,"%s%d_%s_%04d.out",s0ff_base_input,r,kappa_sea,conf_to_read);
                   fin=fopen(s0ft_input,"rb+");
                   if(fin == (FILE *)NULL) {
                        fprintf(stderr,"Unable to open file %s for reading\n",s0ft_input);
                        return EXIT_FAILURE;
                   }
	           for(int imass=0; imass<Nmass; imass++)
                   {
                  // printf("im=%d\n",imass);
		   for (int imom=0; imom<Nmom; imom++)
			for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
			fread(&s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext],sizeof(complex),1,fin);
		   }
		   fclose(fin);
    	
   	   }
	   }
	   printf("====================\n\n\n");

	    ////////////////////////////////////
            /////         CLUSTERIZE      //////
	    ////////////////////////////////////
	    printf("====================\n");
            printf("CLUSTERIZE\n");

	   complex   aux_conf[Nconf], aux_clust[Nclust+1];

          //def
          colorcolorspinspin ****s0FT=malloc(sizeof(colorcolorspinspin***)*(Nclust+1));
          for (int iclust=0; iclust<Nclust+1; iclust++) s0FT[iclust]=malloc(sizeof(colorcolorspinspin**)*(Nmass));
          for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) s0FT[iclust][imass]=malloc(sizeof(colorcolorspinspin*)*(Nmom));
	  for (int iclust=0; iclust<Nclust+1; iclust++)  for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0FT[iclust][imass][imom]=malloc(sizeof(colorcolorspinspin)*(Nr));

	  for (int r=0; r<Nr; r++)
	    for(int imass=0; imass<Nmass; imass++)  for (int imom=0; imom<Nmom; imom++)
            {
		for (int ic_ext=0; ic_ext<3; ic_ext++) for (int is_ext=0; is_ext<4; is_ext++) for (int ic_int=0; ic_int<3; ic_int++) for (int is_int=0; is_int<4; is_int++)
		{	
	               for(int iconf=0;iconf<Nconf;iconf++) aux_conf[iconf]=s0FT_r[iconf][imass][imom][r][ic_int][ic_ext][is_int][is_ext];
                       conf2clust(aux_conf, aux_clust,Nconf,Nclust,Nex);
                       for(int iclust=0;iclust<Nclust+1;iclust++) s0FT[iclust][imass][imom][r][ic_int][ic_ext][is_int][is_ext]=aux_clust[iclust];
		}
            }
	   
            printf("====================\n\n\n");




//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@//

	/////////////////////////////////////////////////////
        // FORM FACTOR OF THE FOURIER-TRANSFORMED PROPAGATOR
	//////////////////////////////////////////////////////

        printf("====================\n");
        printf("FORM FACTOR\n");

	//def
	complex *****s0_ff=malloc(sizeof(complex****)*(Nclust+1));
        for (int iclust=0; iclust<Nclust+1; iclust++) s0_ff[iclust]=malloc(sizeof(complex***)*(Nmass));
        for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) s0_ff[iclust][imass]=malloc(sizeof(complex**)*(Nmom));
	 for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) s0_ff[iclust][imass][imom]=malloc(sizeof(complex*)*(Nr));
	for (int iclust=0; iclust<Nclust+1; iclust++) for(int imass=0; imass<Nmass; imass++) for (int imom=0; imom<Nmom; imom++) for (int r=0; r<Nr; r++) s0_ff[iclust][imass][imom][r]=malloc(sizeof(complex)*(3));
	

	 for (int imom=0; imom<Nmom; imom++)
	 {
		        for (int iclust=0; iclust<Nclust+1; iclust++) for (int imass=0; imass<Nmass; imass++) for (int r=0; r<Nr; r++)
			compute_s0_ff(s0FT[iclust][imass][imom][r],s0_ff[iclust][imass][imom][r],SinP[imom],SinP2[imom]);
	  }

      

        printf("====================\n\n\n");
        

  
///***************************************************************************************************************///////
//********************************************************************************************************************///
///***************************************************************************************************************///////

	/////////////////
	// Write output
	////////////////



     complex **aux_allmomenta=malloc(sizeof(complex*)*(Nmom));
     for (int imom=0; imom<Nmom; imom++) aux_allmomenta[imom]=malloc(sizeof(complex)*(Nclust+1));


      FILE *fout=NULL;      
      char s0ff_david[1024];
      const int  Nlist=4;
      int Listmom[]={100,120,132,140}; 

      printf(" WRITING OUTPUT IN %s\n",dir_output);
      for (int sigma=0; sigma<3; sigma++)
      {

        for (int r=0; r<Nr; r++) for (int imass=0; imass<Nmass; imass++)
        {
        // ff S(p)

         sprintf(s0ff_output,"%s_r%d_m%d_sigma%d.out",s0ff_base_output,r,imass,sigma+1);
         for (int imom=0; imom<Nmom; imom++) for (int iclust=0; iclust<Nclust+1; iclust++) aux_allmomenta[imom][iclust]=s0_ff[iclust][imass][imom][r][sigma];
         write_output_allmomenta(s0ff_output,aux_allmomenta,SinP2,Nmom,Nclust);

         sprintf(file,"ff_r%d_m%d_sigma%d",r,imass,sigma+1);
         write_output_momentabymomenta(dir_output,file,aux_allmomenta,SinP2,Nmom,Nclust);

/*	// David Output
	sprintf(s0ff_david,"%s/FF_SOMEP/ff_r%d_m%d_sigma%d.out",dir_output,r,imass,sigma+1);
	fout=fopen(s0ff_david,"w");
        if(fout == (FILE *)NULL) {
                fprintf(stderr,"Unable to open file %s for writing\n",s0ff_david);
                return EXIT_FAILURE;
         }
         for (int ilist=0; ilist<Nlist; ilist++)
         {
			  int imom=Listmom[ilist]-1;
                          for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=aux_allmomenta[imom][iclust];
                          fprintf(fout,"P%d [imom=%d] %8.10f \t %8.10f %8.10f  + I %8.10f %8.10f \n",Listmom[ilist],imom,SinP2[imom],creal(aux_allmomenta[imom][Nclust]),creal(error_jack(aux_clust,Nclust)),cimag(aux_allmomenta[imom][Nclust]),cimag(error_jack(aux_clust,Nclust)));
        }

        fclose(fout);

*/

        }
      }

////// FINISH
	printf("Free memory\n");

	free(s0_ff); 
	
	printf("Calc_FF_TEST.c finished OK ! \n");
	
  return 0;
}
